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NUMERICAL STUDIES OF UNSTEADY TRANSONIC FLOW OVER AN OSCILLATING AIRFOIL 
W J Chyu and S S Davis 

NASA Ames Research Center, Moffett Field, California 94035, U.S.A. 


SUMMARY 

A finite-difference solution to the Navier-Stokes equations combined with a time-varying grid-generation 
technique was used to compute unsteady transonic flow over an oscillating airfoil. These computations were 
compared with experimental data (obtained at Ames Research Center) which form part of the AGARD standard con- 
figuration for aeroelastic analysis A variety of approximations to the full Navier-Stokes equations was used 
to determine the effect of frequency, shock-wave motion, flow separation, and airfoil geometry on unsteady 
pressures and overall air loads Good agreement is shown between experiment and theory with the limiting 
factor being the lack of a reliable turbulence model for high-Reynolds-number, unsteady transonic flows. 


SYMBOLS 


a 

c 



e 

J 

k 

M 

P 

Pr 

Re 


u,v 


speed of sound 
chord length 

specific heat at constant pressure 
pressure coefficient normalized by 

nth complex component of Cp per 
radian of a 

total energy per unit volume normalized 
by P„a£ 

Jacobian of transformation between physi- 
cal and computational coordinates, 

J = U x n y " 5 y n x) 

reduced frequency, uc/EU^ 

Mach number 

pressure, normalized by p„a* 

Prandtl number, pCp/<* 

Reynolds number, p^a^c/p,,, 

Cartesian velocities, normalized by a,*, 
free-stream velocity 


x.y 


a 

a 

Y 

K 


5>n 

P 

li, X 


T xx’ T yy’ T xy 


physical Cartesian coordinates, normalized 
by c 

angle of attack 

amplitude of airfoil oscillation 

ratio of specific heats 

coefficient of thermal conductivity, 
normalized by <* 

computational coordinates in streamwise 
and normal spatial directions 

density normalized by p„ 

first and second viscosity normalized by 

u* 

time variable, t = ta„/c 
skin friction 


Subscripts 


» free-stream value 

m mean value 

Superscripts 

* reference value 


INTRODUCTION 

The proper computation of unsteady, transonic viscous flows around an oscillating airfoil remains an 
outstanding and important problem in fluid dynamics An efficient and complete analytical capability to 
predict the flow would find immediate applications in the treatment of aeroelastic (flutter and buffet) and 
control-response analyses for both fixed- and rotary-wing aircraft. The theoretical analysis of transonic 
flows is complicated by the presence of mixed subsonic and supersonic regions within the flow field For the 
unsteady flow field, such as that surrounding an oscillating airfoil, additional considerations are needed 
to treat time-dependent aerodynamic loads caused by moving shock waves and boundary-layer interactions At 
Ames Research Center, studies related to these problems have been conducted in recent years, both theoretically 
and experimentally, to predict and clarify many aspects of these flows. 

The physics of unsteady, transonic flow can be simulated at various levels of inviscid and viscous 
approximations For those cases in which viscous effects dominate, computations based on the Navier-Stokes 
equations are needed. Beam and Warming (Ref 1) report on an efficient implicit numerical algorithm for com- 
pressible viscous flow in which an implicit factorization scheme is used Steger (Ref 2) applied this tech- 
nique to the unsteady, compressible Navier-Stokes equations, using the thin-layer approximation In Ref 3 he 

numerically demonstrated the dynamic phenomena of transonic buffet and aileron buzz Recently, Chyu et al 
(Refs 4-7) applied a related numerical method for analyzing unsteady, transonic flows over an oscillating 
airfoil, the method combined the numerical technique reported by Steger with a new and efficient time-varying, 
grid-generation technique suitable for the treatment of moving airfoils. The method was recently applied by 
Horiuti et al (Ref. 8) to the analysis of transonic flows over an airfoil with oscillating flap 

The purpose of this paper is to summarize the results of a 6-year effort in the Ames Aerodynamics Research 
Branch to measure and calculate unsteady transonic flows By using certain data sets identified by AGARD as a 
standard configuration (Ref 9), a series of increasingly complicated unsteady transonic flow cases will be 
analyzed The effects of frequency, flow separation, and airfoil geometry will be studied, and it will be 
shown that increasingly complex flows demand increasingly sophisticated equations to correctly model the flow 
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pattern. The computed results presented here also show that the differences between the full and thin-layer 
Navier-Stokes equations are not significant for this class of flows. 


NUMERICAL PROCEDURE 

A numerical procedure is briefly described for computing the unsteady flow field induced by an oscillating 
airfoil in transonic flow The development of the equations, finite-difference algorithm, boundary treatment, 
and grid system follows closely that developed in Refs. 1, 2, and 6. 

Governing Equations 

An implicit finite-difference scheme was used to solve the full time-dependent compressible Navier-Stokes 
equations in conservation-law form. The numerical algorithm adopted for this study is the Beam-Warming 
spatially factored scheme. The application of the algorithm to the Navier-Stokes equations subject to the 
general transformation, £ = e(x,y,t), n = n(x,y,t), t = t (Fig 1) is present in the Appendix 

For turbulent flows, the viscosity coefficient is computed using the two-layer, Cebeci-type eddy-viscosity 
model reported by Baldwin and Lomax (Ref 10). The main features of the Baldwin-Lomax model are the determina- 
tion of the eddy mixing-length scale based on the local vorticity In the present work, the instantaneous 
vorticity is calculated after each time-step of computation 

Boundary Conditions 

On the outer boundary of the grids, free-stream conditions were applied On the airfoil surface, the 
Cartesian velocities are 


(") - r'P 

-vipo 

\v/ [_' n X 

5«JL» - "tJ 

where the contravariant velocity components U and V 

are defined as 

u - + y + y ■ 

V = n t + n x u 

and 


h = ’Vx ' Vy ' 

n t = ‘Vx ‘ 


with x T and y T determined from local airfoil surface velocities For viscous flow, the no-slip condition 
requires that U = V = 0 at the airfoil surface 

For inviscid computations, a tangency condition is imposed on the airfoil surface by setting V to zero, 
and U is determined at the body surface by a linear extrapolation from the flow field determined at the 
previous step. At the trailing edge of the airfoil, the Cartesian velocities u and v are set equal to x T 
and y T , respectively, to satisfy the Kutta condition. The surface pressure is determined from the normal 
momentum equation in both the inviscid and viscous computations 

Grid Generation 

Stationary grids were first constructed for the airfoil at its extreme angle-of-attack positions Grids 
at intermediate angles of attack were obtained from those at the extreme positions by spatial interpolation 
The extreme-position grids were generated by numerical solutions of elliptic equations (Refs 4 and 11) This 
technique permits grid points to be specified along the entire boundary of the computational plane In Fig 1, 
the boundary is indicated by a curve a-b-c-d-e-f-g-h-i-a in a physical plane that encompasses the airfoil, 
wake, and far-field Application of the method then generates a smoothly spaced, non-overlapping grid at the 
interior points. 

In the present work, boundary points were specified at fixed locations along the airfoil, with grid spac- 
ing clustered near the leading and trailing edges of the airfoil and near the shock-wave position. The grid 
points on the wake (a line, ab or de, in Fig 1) were chosen to lie on a third-degree polynomial arc tangent to 
the bisector of airfoil at its trailing edge and returning to the airfoil centerline at the downstream outer 
boundary This procedure aligns the grid points in the wake with the approximate initial and final directions 
of the wake flow (Fig 2), which was found to improve the efficiency of the flow-field computations. Finally, 
grid points at the free-stream boundary are chosen to lie approximately 8 chord lengths above and below the 
airfoil and 8 chords upstream and downstream of the airfoil leading edge. The outer boundary points remain 
fixed in space for all angles of attack This general approach of fixing the outer boundary points permits 
treating an airfoil oscillating in the proximity of a wall, a second airfoil, or a flap oscillating behind a 
fixed airfoil, as was done in Ref 8. 

Once the grid is specified on its boundaries, an elliptic solver is used to generate a smoothly spaced 
grid at the interior points This grid is then re-spaced, or clustered, along j; = const lines (lines moving 
away from the airfoil), using the weighted coordinate-stretching technique advocated in Ref. 4 

Allowing the airfoil surface to vary within a stationary outer boundary requires that a new grid be 
generated at each time-step of the computation To reduce the computational effort needed to repeatedly 
generate the grid, a novel grid interpolation scheme was devised (Ref. 4). Grids were generated at the 
extreme angle-of-attack positions of the prescribed airfoil motion, using the elliptic equation and the weighted 
coordinate-stretching techniques described previously, these grids were then stored. Grids needed at inter- 
mediate airfoil patterns were found from interpolation along the circular arcs that were assumed to represent 
the locus of the grid-point movement The radius of curvature for each point of the grid was taken as the 
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distance between the grid point and the fixed pivot point about which the airfoil was oscillating. With the 
radius of curvature and the extreme positions known, the center of curvature was then computed for each grid 
point. Grid points for intermediate values of the angle of attack were found from a linear interpolation of 
arc length along the curves defined in this manner. This method of grid interpolation eliminated entanglement 
of coordinate lines during the airfoil motion, even within the dense grids in the boundary-layer region. An 
example of the grid patterns generated by this interpolation method is shown in Fig. 2. Grids similar to those 
used in computations obtained at the extreme angle-of-attack positions are shown in Figs. 2(a) and 2(b), and 
an interpolated grid is shown in Fig 2(c). 

Complex Representation of Surface-Pressure Variation 

The time variation of the surface-pressure distributions is expressed in terms of its Fourier components. 
The harmonically varying angle of attack can be expressed as 

a(t) = a m + Re(o e lut ) 

where Re represents the real part of the complex argument. The Fourier series representation of the surface- 
pressure coefficient is 


C.(x/c,t) = C (x/c) + Y Re[c" -(x/c)5 e lnut ] 
n=i 

where C Pm (x/c) is the mean value of the local surface-pressure coefficient and [C p> 5 (x/c)] is the nth com- 
plex component of the local unsteady pressure coefficient, per radian The real and imaginary value of C p> s 
can be expressed as 

Re[Cp-(x/c)] = (Cp/a)cos * n 
and 

Im[Cp, a ( x /c)] = (Cp/5)sm $ n 

where C p is the nth harmonic unsteady pressure (real) and $ n is the nth harmonic phase shift between 
the angle of attack and the pressure response. 


RESULTS 

Experimental Measurements 

A series of experimental measurements has been carried out in the 11- by 11-Foot Transonic Wind Tunnel at 
Ames Research Center. In these experiments, a NACA 64A010 airfoil and an NLR 7301 airfoil underwent small- 
amplitude harmonic oscillations in pitch about various pivot axes located along the airfoil chord (Refs. 9 
and 12-17) The test Mach numbers ranged from 0.5 to 0.8, chord Reynolds numbers from 2 5 x 10 s to 12 5 * 10 s 
were obtained. Unsteady pressure distributions were measured on the airfoil during the course of the oscilla- 
tions, and unsteady aerodynamic forces and moments were obtained subsequently from integration of the pressure 
distributions. The cases of interest for the present study of the NACA 64A010 airfoil had the airfoil oscil- 
lating about a pivot at x/c = 0 25, with an amplitude of oscillation, 5 = 1°, and a mean value of the pitch 
angle a m = 0° and 4° (With the pitching axes fixed, the pitch angle is equivalent to the instantaneous 
angle of attack of the airfoil ) The range of reduced frequency of the oscillation, k, was taken from 0.025 
to 0 20, and the Mach and Reynolds numbers were held fixed at M„ = 0 8 and Re = 12 x 10 6 , respectively 
(Table 1). 

For the NLR 7301 supercritical airfoil, the airfoil was caused to oscillate in pitch about a pivot at 
x/c =04, with an amplitude of oscillation 5 = 0.5°, and a mean value of the pitch angle, a m = 0 37°. The 
range of reduced frequency k was taken from 0.05 to 0.2, and the Mach and Reynolds numbers were held fixed 
at M„ = 0 75 and Re = 11.4 x 10 6 , respectively (Table 1). 

Computations and Code Performance 

Both inviscid and viscous unsteady computations were carried out for the flow conditions previously out- 
lined Viscous computations were based on both the thin-layer and the full Navier-Stokes equations The 
computational grids have a common outer boundary located 8 chord lengths above and below the airfoil leading 
edge with the airfoil at a = 0. The domain of the grid in the upstream and downstream directions was set to 
8 chord lengths from the airfoil leading edge. 

A grid having 139 x 49 points in the 5 and n directions, respectively, was used throughout the computa- 
tion on the Cray-XMP computer at Ames. The airfoil configurations were described with 117 grid points The 
computations were started from a steady-flow solution obtained at neutral position of airfoil oscillation, 
computations were continued through three cycles of oscillation to reach a periodic solution. 

Including the viscous terms requires additional computational effort per grid point, and the fine mesh, 
needed to resolve the boundary-layer region, necessitated a smaller time-step than the one used in the inviscid 
computations. Decreasing the frequency of the airfoil oscillations also requires an increased number of 
iterations per cycle to maintain the time-accuracy requirement of the Courant number For airfoil oscillations 
at a reduced frequency, k = 0.2, the computation required 2620 iterations per cycle, whereas at k = 0 05 
and 0 025, the computation required 4320 iterations per cycle. With vectorized codes on the Cray-XMP Computer, 
the required computational time per iteration for 139 x 49 grid points was 0.3 sec for the full Navier-Stokes 
code, 0.22 sec for the thin-layer Navier-Stokes code, and 0.17 sec for the inviscid computations. 
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instantaneous Surface-Pressure and Shock-Wave Excursion 


Case NACA 64A010 at low incidence, a = 0° + 1° cos ut, k = 0.025 and 0.20 (Table 1, Fig. 3). 

Computed and experimental surface-pressure distributions, obtained as the airfoil angle of attack varied 
harmonically from 1° to -1° at k = 0 2, are shown in Fig. 3. Computations are shown for the thin-layer form 
of the Navier-Stokes equations. Results are presented for only one-half cycle of oscillation, since the 
motion is symmetrical . At the extreme angle of attack, a = 1°, a supersonic region is evident on the upper 
surface that terminates in a shock at x/c = 0 55. As the angle of attack decreases, the flow speed gradually 
decreases on the upper surface, with corresponding increases in pressure At the same time, the shock strength 
decreases, as evidenced by the decrease of the pressure rise at the shock By contrast, the flow speed on the 
lower surface of the airfoil increases, resulting in reduced surface pressures and a reversed flow pattern 
The comparisons of the computational results (with free-stream boundary condition) and the experimental mea- 
surements show that the computed pressures ahead of the upper-surface shock are underpredicted. The studies 
reported in Refs 18 and 19 and the computation with the experimentally measured wall-pressure boundary- 
condition of Refs 8 and 20 indicate that this underprediction of the pressure is in part a result of the wind- 
tunnel-wall interference The pressure distributions downstream of the shock where the flow is not separated 
are well predicted 

Computed and experimentally measured time-variations of the local pressure coefficient on the upper sur- 
face of the airfoil for a complete cycle of oscillation at k = 0 025 and 0 20 are shown in Fig. 4. In this 
case, results from all three forms of the equation (inviscid, thin-layer, and full Navier-Stokes) are all in 
good agreement The surface pressures display sinusoidal variations ahead of the shock wave (points A, B, 
and C), and nonlinear variation in the shock region (point D). The pressures also show hysteretic variations 
at the higher frequency, k = 0 2. Significant hysteresis is evident at point D in the shock region where 
the pressure undergoes a jump (indicating the passing-over of the shock wave) during decreasing angle of 
attack and smooth-pressure recovery during increasing angle of attack The computational results shown in 
Fig 4 using the inviscid and viscous approximations are all in good agreement with experimental data that 
indicate that viscous effects are small in these flow regimes 

The loci of the computed and measured shock-wave excursions are depicted in Fig 5, which shows that the 
shock moves linearly with the airfoil motion and travels over approximately 71 of the chord with the midchord 
as neutral position (approximately) The trend of the excursions is well predicted by either the inviscid- or 
viscous-flow computations. The results of the viscous computation, however, agree better with the experimental 
data 


Case NACA 64A010 at high incidence, a = 4° + 1° cos wt, k = 0 20 (Table 1, Fig 6). 

This case is different from the previous case of low-incidence flow mainly in respect to the mean angle 
of attack, it is concerned with a flow strongly governed by shock-wave/boundary-layer interactions and a moving 
shock wave with greater strength variation (cf Refs. 16 and 21 for details). The test also showed that for 
k = 0 2, shock-induced separated flow is present for most angles of attack 

Computed surface-pressure distributions based on both of the viscous (full and thin-layer) assumptions are 
shown in Fig 6 as the airfoil angle of attack is varied from 5° to 3° and back to 5° For inviscid computa- 
tions, only typical instantaneous pressures at a = 5° are shown in Fig 6(a), for the pressure distributions 
retain about the same features during the angle-of-attack variation and are not in good agreement with the 
experimental data Also shown in Fig 6 are the experimentally measured pressure distributions. Within this 
angle-of-attack range, a supersonic region is evident only on the upper airfoil surface, whereas, on the lower 
surface the flow remains subsonic. Only a typical lower-surface-pressure distribution is shown in Fig 6(a), 
for it does not vary significantly with airfoil motion 

An examination of the experimental data in Fig 6, starting with the airfoil leading edge on the upper 
surface, shows that the pressure decreases rapidly within a short distance from the leading edge 
(x/c =0-0 05) to a fairly level plateau The computed results of this plateau pressure level by both the 
inviscid and viscous models are all in fairly good agreement with one another, in spite of the large discrep- 
ancy among the theories in the region downstream of the shock wave This indicates that downstream effects, 
including shock-wave discontinuities and shock-wave/boundary-layer interactions have little effect on the 
plateau pressure (or supersonic) region of the airfoil In a manner similar to that of the low-incidence case, 
the computed pressures ahead of the shock are underpredicted in part because of the wind-tunnel-wall inter- 
ference Recent experiments show that the interference was more significant as the incidence of the flow to 
the airfoil is increased (Refs 18 and 19). 

The comparisons of the inviscid computations and the experimental measurements show that the shock 
strength is overestimated and that the computed shock wave is positioned too far downstream when the inviscid 
theory is used The computed results with the full and the thin-layer Navier-Stokes equations showed no sig- 
nificant difference in the magnitude of instantaneous surface pressures Although the location of the shock 
wave differs by about 0%-3% of the chord, neither computation gave consistently closer agreement with experi- 
mentally measured surface pressures over the entire range of the airfoil oscillation. 

Also shown in Fig 6 is the measured upper-surface-pressure distribution, it shows an aft-shock pressure 
recovery A slow pressure recovery is typified by a pressure distribution with a large bump behind the shock 
wave, such as those shown in Fig 6(a) for a = 5°, a fast recovery is typified by a smooth pressure increase, 
such as those for a = 3 13° (Fig. 6(h)). The aft-shock pressure bump is indicative of probable flow separa- 
tion induced by the shock-wave/boundary-layer interactions, whereas the fast pressure recovery indicates 
probable attached flow on the airfoil surface behind the shock wave. Steady flow interferograms obtained by 
Johnson and Bachalo (Ref 21) on the same airfoil under the same flow condition show extensive shock-induced 
separation for a = 5° and attached flow for a - 3° (unpublished data) in the aft shock region. The measured 
pressure variations aft of the shock wave in Fig 6 show that a slow pressure recovery (indicative of separated 
flow) is maintained for most of the angle-of-attack variation, except in the upward motion of the airfoil when 
the angle of attack is varied from 3° to 3 13° Fairly good agreement between the computed and the measured 
pressure recovery was obtained with the Navier-Stokes equations Although the transition angle of attack 
between the slow and fast pressure recoveries (indicative of separated and attached flow, respectively) is not 
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accurately predicted, some of the Important effects of the shock-induced separation and the shock-wave/ 
boundary-layer interactions are clearly demonstrated by the aft-shock pressure distribution in Fig 6. 
Typically, Fig. 6(b) depicts a slow aft-shock pressure recovery at a = 4 87°. The corresponding computed 
instantaneous velocity vectors depicted in Fig. 7 show shock-induced separation at the same angle of 
attack. 

Figure 8 shows Instantaneous pressure distributions on the airfoil upper surface when the angle of attack 
is increasing or decreasing and passing through the neutral position a = 4°. The results exhibit a substan- 
tial phase shift of the pressure relative to the airfoil motion, the shock positions differ by about 10% of the 
chord, and the flow aft of the shock is separated during decreasing angle-of-attack motion and is attached 
during increasing angle-of-attack motion. 

The loci of the computed and measured shock-wave excursions are depicted in Fig. 9, which shows that the 
shock travels over approximately 13% of the chord length, with the neutral position at about x/c = 0.45. 

In this figure, some of the important features of unsteady transonic flow at high incidence are depicted. For 
a flow without separation, such as the case of low flow-incidence (Fig. 5), and in accordance with inviscid 
theories, the shock wave moves downstream with increasing angle of attack. In this case, however, the shock 
wave initially moves downstream with increasing angle of attack (from 3° to 5°), but starts to move upstream 
at 4.8° owing to the occurrence of flow separation, as was indicated in Fig. 6(1) and 6(a) by the appearance 
of slow aft-shock pressure recovery. It appears that this important "retrograde" shock motion cannot be pre- 
dicted with inviscid transonic theories. 

The results of the "thin-layer" computations shown in Ref. 6 differ somewhat from those shown in Fig. 6. 

In the present computation, grid spacing was densely clustered near the airfoil surface (about two grid points 
in the sublayer) to more accurately compute the eddy-viscosity terms in the turbulence model. 

Case NLR 7301 at low incidence, a = 0 37° + 0.5° cos ait, k = 0.05 and 0.20 (Table 1, Fig. 10). 

Figure 10 shows data for a flow that differs from the previously considered flow in terms of airfoil 
geometry, pitch-axis location, and amplitude of oscillation. The NLR 7301 airfoil was designed for "shock- 
free" flow and is an example of a relatively thick (t/c = 0.165) modern supercritical airfoil (Refs. 22 
and 23) Experimental data show that the flow was shock-free at = 0.75 as the airfoil was held fixed at 
a = 0.37°. At off-design condition, as in this case, experimental data in Fig. 10 show that a weak shock was 
evident on the airfoil upper surface. 

Both the computed and the experimentally measured pressures show small variations, relative to the cases 
of the NACA 64A010 airfoil previously studied, because of the oscillatory motion. Thus, only the instantaneous 
pressures at the extremes of angle of attack are shown in Fig. 10. The pressure variations in Fig. 10 show a 
pressure rise to a plateau level immediately after the leading edge owing to the blunt and thick leading sec- 
tion of the supercritical airfoil. Unlike the conventional airfoil, wavy pressure distributions are shown in 
Fig 10 in the plateau region (x/c = 0.1 - 0 6). The time-variation of the computed and measured local pres- 
sures in the plateau region also show small unsteady pressure fluctuations (or higher harmonic component) in 
addition to the overall pressure variation, indicating that the supercritical airfoil is very sensitive to the 
unsteady effects from the flow field. Although these small pressure variations in the plateau region are not 
computationally well predicted, the computational results in Fig. 10 show the general feature of the pressure 
typical to the supercritical airfoil. Comparison of pressures at the extremes of angle of attack between the 
two frequencies (k = 0.5 and 0 20) indicates that the pressures are less sensitive to the airfoil motion at 
higher frequency. At lower frequency, k = 0.05, only those pressures in the plateau region are affected by the 
airfoil motion. 

Harmonic Analysis of Unsteady Pressures 

Both computed and measured unsteady pressures have been expressed in Fourier components up to and includ- 
ing the third harmonic The higher harmonic components are not presented in the following figures since the 
modal content decreases rapidly with increasing mode number Distinctive features of the pressures depicted in 
Figs. 11-14 are discussed in the following paragraphs 

Case. NACA 64A010 at low incidence, a = 0° + 1° cos u>t, k = 0.025 and 0.20 (Table 1, Fig. 11) 

Illustrated here are the basic characteristics of the unsteady surface pressure for the conventional 
airfoil as it is subjected to a flow with small flow-incidence or a weak shock-wave/boundary-layer interaction 
or both. The mean-pressure features a distribution similar to those of steady flow, and is reasonably well 
predicted except in the shock region, where the viscous computations show closer agreement with the experimen- 
tal data. 

Examination of the higher frequency (k = 0.2), in-phase (Re) and out-of-phase (Im) components indicates 
that the pressures forward of the shock wave (x/c < 0 5) contain both components in about equal magnitude. 
Pressures behind the shock wave, however, contain mostly out-of-phase components. In the low-frequency case 
(k = 0.025), pressures in front of the shock wave contain more in-phase than out-of-phase contributions, 
whereas behind the shock wave both components are small and of about the same magnitude. The out-of-phase 
components for both frequencies, k = 0.025 and 0.2, also show that the component in front of the shock is 180° 
out-of-phase with those behind the shock, and a rapid phase shift abruptly takes place at the shock. At both 
frequencies, the real and imaginary components shown in Fig. „ 11 and the measured data for other frequencies 
(Ref. 9) exhibit similar distributions over the airfoil surface and indicate that the variation of the compo- 
nents is approximately linear with the frequency in the flow for which shock-wave/boundary-layer interactions 
are minimal . 

Case. NACA 64A010 at high incidence, a = 4° + 1° cos u>t, k = 0.05 and 0.20 (Table 1, Fig. 12). 

The high-incidence flow (am = 4°) considered here differs from the low-incidence case (a m = 0°) in the 

mean angle of attack and is characterized by dominant shock-wave/boundary-layer interactions. The mean 
unsteady pressure shown in Fig. 12 exhibits a rapid decrease from the airfoil leading edge to a plateau pres- 
sure level ahead of the shock. It also shows an abrupt pressure jump at the shock region followed by a slow 
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pressure recovery (indicative of separated flow). Only the mean values of the pressure computed with inviscid 
assumption for k = 0.2 are depicted in Fig 12. It shows that the predicted shock position is too far 
downstream. 

Examination of the harmonic components of the pressure for the high-frequency case, k = 0.2 in Fig. 12, 
indicates that the upper-surface pressure in front of the shock wave contains both in-phase (real part) and 
out-of-phase (imaginary part) components of about equal magnitude, however, the pressure immediately behind 
the shock wave (x/c = 0.6) contains mostly in-phase components Toward the airfoil trailing edge 
(x/c = 0.6 - 0.9), the pressure increases gradually in its in-phase and rapidly in its out-of-phase components 
In the trailing-edge region, the pressures contain mostly out-of-phase components. The variation to these 
harmonic components behind the shock distinguishes the present case from the previous case (low-incidence flow) 
'here the flow behind the shock wave is attached and the pressure contains only the out-of-phase components 

The effects of frequency can be evaluated by comparing the harmonic components of the pressure for both 
k = 0 05 and 0 2 in Fig. 12. This comparison indicates that the surface pressure at the lower frequency 
k = 0 05 contains mostly in-phase components. In the shock region, the in-phase components at the lower 
frequency are 180° out of phase with those at the higher frequency Toward the airfoil trailing edge, the 
'n-phase components rapidly decrease at the lower frequency, whereas they exhibit the opposite trend (gradu- 
ally decrease) at the higher frequency. These nonlinear variations of the harmonic components with frequency 
differ distinctly from those previously considered in low-incidence flow, particularly in the shock wave and in 
the aft-shock separated regions of airfoil surface They also show strong nonlinear variations with frequency, 
thus reflecting the complex effects from shock-wave/boundary-layer interactions and shock-induced separation 

Pressure variations on the lower airfoil surface in Fig. 12 do not exhibit variations as rapid as those 
on the upper surface. The harmonic components of the pressure show a nearly linear variation from the leading 
t- the trailing edges of the airfoil 

Case NLR 7301 at low incidence, a = 0 37° + 0.5° cos wt, k = 0 05 and 0.20 (Table 1, Figs. 13-14). 

The case considered here differs from the previous studies in the airfoil geometry Figure 13 shows a 
plateau in the mean pressure that is reached rapidly from the blunt leading edge of the thick airfoil and that 
is sustained up to the region of rapid compression (but not a shock) on the airfoil (x/c = 0 6). The mean 
unsteady pressure in the aft-region of compression exhibits a rapid pressure recovery that indicates an unsepa- 
rated flow 

The harmonic components of pressure show that the pressure in the fore section of the airfoil contains 
both in-phase and out-of-phase components of about the same magnitude (as in the NACA 64A010 in high incidence 
flow, and at k = 0 2 only. Fig 12) The magnitude of these components is significantly increased at about 
60 % of the chord, a result of the occurrence of rapid compression. 

Comparison of the harmonic components of the pressure for both frequencies, k = 0 2 and 0 05, indicates 
that the out-of-phase components do not diminish at a lower frequency, but remain at about the same magnitude 
in both frequencies This is contrary to the case of the conventional airfoil (Figs 11 and 12). The in-phase 
components, however, show an increase in magnitude at a lower frequency that is comparable to that seen in the 
case of the conventional airfoil. 

The computed magnitudes of the harmonic components (Fig 13) exhibit only the qualitative trends of the 
experimentally determined pressure variations However, it should be noted that the computation assumed 
uniform free-stream boundary conditions and the wind-tunnel -wal 1 effects were not taken into account Measured 
pressure signatures, on the other hand, can be affected by the wind-tunnel-wall interference (Ref 8) and by a 
small change in free-stream turbulence level in different wind tunnels (Refs 16 and 22) Experiments show 
that the pressure is also sensitive to geometry even at such small protuberances as a pressure transducer in 
the shock region The mathematically smoothed, measured airfoil profile used in the present computation could 
lead to the discrepancies. 

The phases of the complex pressure components for both frequencies, k = 0 2 and 0 05, are shown in 
Fig. 14, good agreement with experimental data is shown The figure also shows that phase varied gradually in 
the leading section of the airfoil and jumped abruptly at x/c =06. 

Surface pressures on the airfoil lower surface were not experimentally measured However, the computa- 
tional results shown in Fig. 13 show that the harmonic components are minimal on the rear half (concaved por- 
tion) of the airfoil lower surface for the flow conditions investigated. 

The overall pressure distributions displayed in Figs. 13-14 show that the unsteady air loads on the 
supercritical airfoil are distinctly different from those on the conventional airfoil 


CONCLUSION 


The ability of a Navier-Stokes code to compute unsteady transonic flow over an airfoil in oscillatory 
motion has been investigated in conjunction with a series of tests conducted in the 11- by 11-Foot Transonic 
Wind Tunnel at Ames Research Center. Two airfoils were considered, a conventional (NACA 64A010) and a super- 
critical (NLR 7301) section. The study of the conventional airfoil disclosed that the unsteady flow field 
and the related aerodynamic load, both in cases of low and moderately high flow incidence, were fairly well 
predicted The computed results successfully and distinctly demonstrated the nonlinear aerodynamics character- 
ized by the shock-wave/boundary-layer interactions and frequency of the airfoil motion The inviscid Euler 
code was shown to be adequate only in the case of low-incidence flow, or for the flow in which viscous effects 
are minimal For the supercritical airfoil, only the computed mean unsteady pressures and phase of the 
Fourier components are in good agreement with the experimental data The harmonic magnitudes of the pressure, 
which are experimentally, as well as theoretically, shown to be extremely sensitive in the transonic regime 
for supercritical airfoils, are predicted in a qualitative manner only The difficulties in correlating com- 

puted results with experimental data in this case could be attributed to the sensitive characteristics of the 
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supercritical airfoil. In summary, the Navler-Stokes code demonstrated fairly well its capability of modeling 

the nonlinear, unsteady aerodynamics for a supercritical airfoil. 

The Navier-Stokes computations with and without the thin-layer approximation showed no noticeable differ- 
ence in the magnitude of instantaneous airfoil surface pressures, except in the location of the shock wave, 
which differs by Q%-3% of the chord. Neither computation gave consistently better or closer agreement with 
experimentally measured surface pressure. A substantial difference in shock position was shown in Ref. 6, 
where insufficient clustering of the grid was made in the sublayer of the turbulent boundary layer in the 
thin-layer Navier-Stokes computations. 

The present code uses an eddy-viscosity model for the turbulent boundary layer that was developed using 
simple steady-flow experiments. For unsteady, high-incidence flow for which viscous effects dominate, an 
improved turbulence model is surely required. However, the demonstrated capability of the code has important 
implications for applications in aeroelastic and control -response analyses and in the study of wake-airfoil 
interactions. 
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APPENDIX 


The Navier-Stokes equations in Cartesian coordinates can be written as (Ref. 2) 

3 t q + 3 X E + a y F = Re _1 (3 x R + 3 y S) (Al) 

where 



where 


T xx = (X + 2u)u x + xv y 

T xy = “< u y + v x> 
t yy = (X + 2y)v y + Xu x 

R “ = UT xx + VT X y + KPr ' 1 ^ ' l)' l 9 x a * 

S " = UT xy + VT yy + lcPr ' 1 ('T * D' l 3 y a 2 
P = (y - l)[e - | p(u 2 + v 2 )J 
a 2 = y(y - l)[f - \ (u 2 + v 2 )] 

and the Stokes hypothesis X + (2/3)y = 0 is assumed. Equation (Al), subject to the general transformation, 
C * C(x,y,t) , n = n(x,y,t) , and t = t, gives 

3 t q + 3 ? (E - R) + 3 n (F - S) = 0 (A2) 


where 


q = q/J 

£ = (C t q + C X E + £ y F)/J 
F - (n t q + n x E + n y F)/J 

and chain rule is applied to stress terms such as x xx * (x + 2y)(c x Ur + n x u n ) + X(c y v^ + ny v n ) , etc. The 
viscous dissipation vectors in Eq (A2) are ' 

R = J‘ l Re* 1 (c x R + c y S) 

§ = J' 1 Re" 1 (n x R + n y S) 

where 

R = ^ 0 ,T xx’ x xy ,R ^ T * ^ = tO»t xy .T yy »S 4 ] T 

The viscous vectors, R and S, contain terms that are functions of (q,q^) and (q.q^), and are written as 

R(q.q^.q n ) = R^q.q^) + ft 2 (q ,q n ) 


where 


S(q,q 5 > q n ) = S^q.q^) + § 2 (q,q n ) 

0 

a i u 5 + “3 V (r 
a 3 U 5 + a 2 V ? 

f (u 2 ) 5 + ^ (v 2 ) 5 + a 3 (uv) 5 + a„( Y - l)- 1 (a 2 ) e 


6 2 (q.q 5 ) = Re’ 1 J' 1 
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R 2 (q.q n ) = Re' 1 J' 1 


^(q.qj = Re" 1 J’ 1 


Yi u n + Y 3 % 

Y - U n + Y * v „ 

T (u2) n + T" (v2) n + y > u \ + Y " vu n + y s (y * 


Y i U C + Y - V C 
y 3 u 5 + V 5 


Y 1 


T + T ( y2 ^ + Y t uv s + y 3 vu s + Y s< Y ‘ 


§ 2 (q.q n ) = Re" 1 J” 1 


s.u + B,v 
in 3 r 


B.u + B,v 
3 n 2 n 


Bi 


y- (u 2 ) n + T (v 2 ) n + B 3 (uv) n + B,(y - D- 1 (a 2 ) n 


Here, 


«i = ^y + K> 

°2 ■ >*(«; + I ty 


a 3 = 7 Vy 

a 4 = KPr'MsJ + c 2 ) 


B i = w(n y + I n x ) 
B 2 = •'(n 2 + j n 2 ) 


B 3 = 3 Vy 


Y i = ^Vy + I 


Y 2 = ^ x n x + 3 S y n y ) 


Y 3 = “(«y n x * 3 Vy* 


S 4 = KPr'Mn 2 + n 2 ) y 4 = uU x n y - § 5 y n x ) 

Y s " ‘^Vx + y y > 

In terms of the viscous terms R x , R 2 , § 2 , S 2 , Eq (A2) is reduced to 
3 T q + S^E + 3 n F = + R 2 ) + 3 n (S 1 + S 2 ) 

A solution to Eq (A3) can be obtained by using a single-step temporal scheme (Ref. 1) 


a ^ n = frr 3 t 


n + rh y" + r^T A ^ n + °[(® - \ - b ) At2 + At3 


(A3) 


(A4) 


where Aq 11 = q n 1 - q , q = q(n At), and n denotes the nth time step of computations The implicit 
trapezoidal scheme takes 0 = 1/2 and B = 0 in (A4) which results in truncation error 0(at 3 ), and the 
implicit Euler scheme takes 0=1 and B = 0 which is 0(At 2 ). The implicit trapezoidal scheme was used lr 
the present computations. 

Applying the scheme (A4) into Eq (A3), one obtains a solution to the Navier-Stokes equations 
Aq" = [3 ^-aE" + AR? + ARj) + 3 n (-AF n + ASj + AS?)] + [^(-E + R x + R 2 )" 

+ 3n(-P + §j + S 2 ) n ] + Aq"* 1 + O[(0 - \ - B)At 2 + At 3 ] (A5) 

where AE n = E n+1 - E n and E n+1 = E(q n+1 ), etc The finite-difference form of q n is q 3>J at x = l ax 
and y = l Ay 

The flux vector increments (aE", AF n , aR 2 , aR 2 , aSj, aS 2 ) are nonlinear functions of q, and are linear- 
ized by using Taylor series expansion. 


AE n = A n Aq" + 0(at 2 ) 


aP" = §" Aq" + 0(AT 2 ) 

where A n = (3E/3q) n , B n = (aF/aq)", or in detail 


(A6) 
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k 0 

1 k i | 

k 2 

1 0 

-u(kjU + k 2 v) + k 2 * 2 

| -(y - 3 )k 2 u + k 0 + k 2 v | 

-(y - l)k 2 v + k 2 u 

1 (y - Dkx 

-v(k 2 u + kjv) + k 2 d > 2 

| -(y - l)k 2 u + k 2 v | 

-(y - 3 )k 2 v + k„ + k 2 u 

I (y - l)k* 

(k 2 u + k 2 v)^- + 2 * 2 ^ 

| 

(? - 

1 y(k 2 u + k 2 v) + k 0 


[ - (y- l)(k 2 u + k 2 v)u , 

- (y - l)(k 2 u + k 2 v)v 

1 


where <j> 2 = (l/2)(y - 1 ) (u 2 + v 2 ), and 



Likewise, the linealization for aR 2 and aS 2 gives 

AR 2 = P" Aq" + R" Aq" + 0(at 2 ) 

AS" = Q n Aq" + C n Aq" + 0(At 2 ) 


where p" = (3R / 3q ) n , R" = (3R 2 /3q )", 0" = ( 3§ 2 / aq ) n , and L n » (3S 2 /3q ) n The linearized ar" and as" in 
the above can further reduced to 5 n 


AR" = (P - K^)" Aq" + 3 ? (R Aq)" + 0 (at 2 ) 
A§" - (Q - L n ) n Aq" + 3 n (C Aq)" + 0 (at 2 ) 


(A7) 


where R^ = a^R and L n = 3 n L n . The viscous matrices R, L, P - , and Q - L n are presented below. 



r 0 i 

0 1 

0 

O' 


-hjU - h 3 v | 

hi 1 

h 3 

0 

R, L, (P - R 5 ), or (Q-L n ) = Re'V 1 

l 

\jr 
I c 
1 1 

1 < 

1 

1 

1 

1 

1 

h, 1 

h 2 

0 


hj-f + u 2 + V 2 ) (h 2 

- hJu + h 3 v 1 (h 2 

- hjv + h 3 u 

\ 


- h 2 u 2 - 2h 3 uv - h 2 v 2 | 

1 

1 




where 


L = °1 

• h 2 

- a 2 , 

h 3 = a 3 , 

h - = “i. 

for 

K , 

L = Bx 

> h 2 

= S 2 , 

h 3 = 6 3 , 

\ = B 4 

for 

L , 

3a l 

35 ’ 

h 2 - 

3a 2 
35 ’ 

05 1^° 

II 

-C 

3a 4 

hl * ~ ~dlT 

for 

P - 

3S 2 
3n 1 

- 

36 2 
3n ’ 

=rU 

11 

-C 

h - ^ 

4 3n 

for 

Q - 


The spatial derivatives 3R 2 /3t and 3§j/3n involve cross derivatives 3 2 ( )/353n, in order to use a 

spatially factorized scheme, these terms are time-lagged as 


aR 2 = ar"' 1 + 0(at 2 ) 


AS" = as"' 1 + 0(at 2 ) 


(A8) 


Substitution from the linearized flux vector increments AE n , aP", arJ, as" { A6 and A7) and the time lagged 
flux vector increments AfiJ and AS" (A8) to Eq. (A5) gives 

{1 * [3, (A - P + R/ - 3 2 R" + 3 n (B - Q ♦ L n )" - 3 2 L n ]}Aq" 


“ T^TT [3 t ( ' § + *i + **>" + 3 n ( - ? + § i + §2>n] + TTT C* E (^2) n-1 + 3 n (AS 1 )"' 1 ] 


y!-jj Aq"' 1 + o[(e - j - e^AT 2 ,(e - 6)at 2 ,at 3 J 


(A9) 
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Here, the notation of the form [3 5 (A-P + Kr)] n aq n implies 3r[(A - P + R^) n _aq n ], etc For second order 
accurate scheme e should be set to 6, and for first order accurate scheme e be set to zero. 

The thin-layer approximation neglects terms 3rK n , 3 r (R, + R,), 3..S, , -P + R f , -Q + L n , and 
Mar,)"- 1 ♦ 3,(15,)^ in Eq. (A9) . 5 5 1 2 ” 1 5 ° 


Implementation of a spatially factored scheme in the left hand side of Eq. (A9) gives 

{i + fvT + V" ' 3 ^ n] } + f+T [3 n ( ® ‘ ^ + C n )n ' 3 n‘- r ’ ] }^ n 3 LHS < A9 > + 0 (at 3 ) 

where LHS (A9) denotes the left hand side of Eq. (A9). 


The following computational sequence is 

{' * T-rs IV* 

{■ * m t>„<s 


used to obtain the Navier-Stokes solutions 

- P + R c ) n - 3*K n ]^Aq n = RHS(A9) 

- 0 + C n ) n - 3 ^ n ]Jaq n = 7? 


r 1 


= q + Aq 


(A10) 


Equation (A10) can be simplified by assuming that the transport coefficients are locally constant, or 
-P + K- = -Q + L n = 0 

The simplified form of Eq. (A10) is thus 

(i + y~r [3. (A - 3^R) n ]|iq n = RHS(A9) 

{l - ^< § - 


In the present computation, the time-lagged viscous term 3^(AR 2 ) n ' 1 + S^tASj)"" was neglected, reducing 
the order of the time-accuracy in the computation of dissipation terms from second to first order The final 
numerical algorithm constructed for the computations is thus first-order time-accurate for the dissipation and 
second-order time-accurate for the convection terms in the Navier-Stokes equations. The spatial derivatives 
in the equations were kept second-order accurate. This brief derivation differs from those of previous 
investigations by including all the relevant terms of the Navier-Stokes equations (Eq (Al)) 



TABLE 1. EXPERIMENTAL DATA THAT WILL BE COMPARED WITH COMPUTATIONAL RESULTS 



Fig. 1. Physical and computational planes. 
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NACA 64A010 

= 0 8, Re = 12 x 10 6 , K = 02 
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EXPERIMENT 


THIN LAYER NAVIER STOKES 




-1 

a, deg 


-1 

a, deg 


(a) K = 0 025 


(b) K = 0 20 


Figure 4 Time-histories of upper-surface pressures NACA 64A010, a = 0° + 1° cos wt 





3-16 


NACA 64A010 
UPPER SURFACE 

M oo = 0 8, Re = 12 x 10 6 , K = 0 2 
EULER 

THIN LAYER NAVIER STOKES 



1 1 I I I 

10-101 
a, deg 


Fig 5 Shock-wave locus on upper surface of the airfoil, a = 0° + 1° cos cot. 



(i) a = 4 0 


(j) o = 3 5 


(k) o = 45 


(I) o = 4 87 


Fig 6. Surface-pressure variation with angle of attack, a = 4° + 1° cos wt 


= 0 8 Re = 12 x 10 6 K = 0 2 



Fig. 7. Computed instantaneous velocity profiles a = 4 87° (a decreasing) 
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Fig 9 Shock-wave locus on upper surface of the 
airfoil a = 4° + 1° cos ut. 
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Fig 10. Instantaneous upper-surface 
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Fig. 13. Mean and first harmonic comply 
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